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I. INTRODUCTION 


There are many flare forecasting models. For an excellent 
review and comparison of some of them see Barnes et al. 
(2016). All these models are successful to some degree, but 
there is a need for better models. We claim the most suc- 
cessful models explicitly or implicitly base their forecasts on 
various estimates of components of the photospheric current 
density J, based on observations of the photospheric mag- 
netic field B. However, none of the models we are aware 
of compute the complete J. We seek to develop a better 
model based on computing the complete photospheric J. 
Initial results from this model are presented in this talk. We 
present a data driven, near photospheric, 3 D, non-force free 
magnetohydrodynamic (MHD) model that computes time 
series of the total J, and associated resistive heating rate 
in each pixel at the photosphere in the neutral line regions 
(NLRs) of 14 active regions (ARs). The model is driven by 
time series of B measured by the Helioseismic & Magnetic 
Imager (HMI) on the Solar Dynamics Observatory (SDO) 
satellite. Spurious Doppler periods due to SDO orbital mo- 
tion are filtered out of the time series of B in every AR 
pixel. Errors in B due to these periods can be significant. 
For the NLR integrated resistive heating rate Q, the num- 
ber of occurrences N(q) of values of Q > q is the cumulative 
distribution function (CDF) of Q. For each AR time series 
it is found to be a scale invariant power law distribution, 
N(Q) « Q7*, above an AR dependent threshold value of Q, 
where S varies little between ARs (§V). For coronal flares, 
N(E), where F is the total magnetic energy converted into 
particle energy, is known from observations to have the same 
form: N(E) «x E~*, s varies between ARs, and the ranges 
of S and s strongly overlap (§V). This strong similarity be- 
tween N(Q) and N(E) suggests the same process that pow- 
ers coronal flares also powers the photospheric Q. Model 
results also suggest it is plausible that the times of large 
spikes in Q, several orders of magnitude above background 
values, are correlated with the subsequent occurrence of M 
or X flares (§IV). These spikes typically occur a few hours 
to a few days prior to M or X flares. The spikes correspond 
to large vertical derivatives of the horizontal magnetic field, 
suggesting strong heating in horizontal current sheets. The 


subset of spikes analyzed at the pixel level are found to oc- 
cur on HMI and granulation scales of 1 arcsec and 12 min- 
utes, suggesting the current sheets are granulation scale. 
Spikes are found in ARs with and without M or X flares, 
and outside as well as inside NLRs, but the largest spikes 
are localized in the NLRs of ARs with M or X flares. Time 
series from more ARs must be analyzed to statistically de- 
termine if this suggested correlation between Q and flares 
is real. 


Il. THE MODEL 


The model is an MHD model that uses all of the available 
HMI B data for a given AR. The model is valid close to 
the photosphere in that the analytical, semi-empirically 
determined expression for B is valid through order z?, and 
the V-B = 0 condition is satisfied through order z. The 
only differential equation in the model is V x A = B, for 
the vector potential A, with the gauge condition V- A = 0. 
This equation is solved analytically. The electric field 
E ~ —c7!0A/0t. The Ohm’s law is E+ (V x B)/c = nJ. 
V is the velocity of the plasma. 7 = 2 x 10~™ s is the 
resistivity. There are no equations for density, temperature, 
pressure, or momentum. Cartesian coordinates (x, y, z) are 
used, with z the height above the photosphere. A data 
determined length scale L(x,y,t) determines 0B/0z at 
z = 0, allowing the complete J(x,y,0,t) to be computed. 
It is enforcing the V-B = 0 condition that allows D(z, y, t), 
and hence J(x, y,0,t) to be computed. The HMI pixel side 
length is A = 0.5”. Pixel overlap results in a magnetic field 
resolution ~ 1’"(~ 725 km). HMI provides time averaged 
data with a resolution of 12 minutes. Data from the 
hmi.sharp_720s_cea data series are used to minimize effects 
of noise. Every 12 minutes HMI provides a full disk map 
of B. HMI B is a function of x,y, and t. Using the CEA 
(cylindrical equal area) data helps minimize projection 
error (Sun 2013). 


A. Magnetic Field 


Let (Lz, L,) be the (x, y) dimensions of the rectangular 
region used to enclose the AR modeled. The number of HMI 


data points covering this region is N = (Nz + 1)(N, + 1), 
where N, = L,/A,N, = L,/A, and Nz, Ny are given by 
the HMI datasets. For any function f(x, y, z,t), define fi. = 
Of /Ox, and similarly for derivatives with respect to y, z, and 
t. Let 


B(x, y, z,t) = e72/L(2,y,t) 3 > Drm(t e2ri( (pe+ +eH) 
n=0 m=0 
(1) 
Here the Dbnm(t) are complex, and L(z,y,t) = 


Lo(a,y,t) + z£[y1(x,y,t)/Lo where Lo and Ly are real 
and determined by the HMI data and the V-B = 0 
condition. Equation (1) is valid for sufficiently small z. 
For z = 0, and given the N vectors B(2;,y;,0,t;) from 
the HMI data for each j, Eq. (1) represents N equations 
for each component of by ,(t). The time series of by (t)is 
determined by performing an FFT on Eq. (1). The 
imaginary part of B must be zero, which is used as a check 
on the numerical solution for the Drm (t). 


B. The V-B=0 Condition 
Define Bo = B(z, y, 0,t). Take the divergence of the real 


part of Eq. (1) and set it equal to zero. Solving the resulting 
equation through order z gives 


Boz 
LN) =a (2) 
x UY 
and 
L 
Ii(2,y,t) = a. (Box Lo,x + BoyLo,y) : (3) 


The right hand sides of Eqs. (2) and (3) are evaluated at 
z = 0. Therefore, Lo and L, are completely determined 
by the HMI data. The resulting expression for B is valid 
through order 27. It is the V-B = 0 condition plus the 
HMI data that determine the z dependence of the model, 
which is what allows the complete expressions for J, and 
Jy to be computed. Without a determination of the height 
dependence of B, the components of J; and Jy that involve 
By, , and Bz, at z = 0 cannot be computed. 


Ill. REMOVAL OF SPURIOUS DOPPLER PERIODS 
FROM THE HMI B DATA 


There is spurious, Doppler shift generated noise in the 
form of 6, 12, and 24 hour period oscillations in the com- 
ponents of B for each pixel, corresponding to frequencies of 
(4.6296, 2.3148, 1.1574) x 10~° Hz. This noise causes a slow 
change in B relative to the granule turnover time since the 
oscillation periods correspond to 20-60 turnover times. The 
model removes this noise from the time series of HMI B for 
each pixel using an FFT based bandpass filter (S14) that 
removes frequencies in an interval of length 0.4166 x 107° 
Hz centered on each of the three noise frequencies.! 


1 An alternative method for mitigating the spurious Doppler effects 
is presented in Schuck et al. (2016) 
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Figure 1. Comparison of Filtered and Un-filtered B, and B, in 
a pixel. By, not shown, is similar to Bz. 


Denote the filtered and un-filtered time series for B, as 
Bzf and Bz,, and similarly for other quantities. Figure 
1 shows the filtered and un-filtered HMI time series of B, 
and B,, the difference between the un-filtered and filtered 
time series, and the magnitude of their ratio for a randomly 
selected pixel from the NLR of NOAA AR 1166. This AR 
is one of the strongly flaring (SF) ARs analyzed here. An 
SF AR is one with M or X flares. A control AR (C AR) 
is one with lower class or no flares. The figure shows that 
the difference between the filtered and un-filtered time se- 
ries of B is significant, so Doppler noise can be significant 
at the single pixel level. Doppler noise can also be signifi- 
cant in quantities that are integrals of pixel level quantities 
over NLRs. For example, again consider the time series for 
NOAA AR 1166 used for Fig. 1. Figure 2 shows the results 
of integrating the filtered and un-filtered pixel level results 
for Q = nJ? and B?/87 over the NLR at each time. The 70 
hour long time interval includes 1 X, 2 M, and 9 C flares. 
For these and subsequent figures, the red, green, and light 
blue vertical lines and their labels indicate the times and 
magnitudes of X, M, and C flares. During the 70 hour time 
interval the number of pixels in the NLR varies across the 
range of ~ 3—6 x 10+. The figures show that the sum of 
un-filtered quantities from each pixel causes a large error in 
the result. 
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Figure 2. Comparison of filtered and un-filtered NLR integrated 
Q=nJ? and B?/8xr. 


IV. HEATING SPIKES IN SF ARs 


In this and the remaining sections of this abstract, unless 
otherwise indicated, Q stands for the spatial integral of the 
pixel level Q over the NLR. Figures 3 and 4 show the time 
series of Q for the 7SF ARs. The times of the largest spikes 
in Q relative to the times of M and X flares in these plots 
indicate it is plausible that these spikes are correlated with 
the subsequent occurrence of M or X flares, and indicate the 
need to analyze time series of Q for more ARs to statistically 
determine if such a correlation exists. The plausibility of 
a correlation is suggested by the following trends shown in 
the plots. For AR 1158 the largest spike by a factor ~ 25 
occurs ~ 38 — 68 hours before the X and M flares. For AR 
1166 the largest spike is ~ 3.5 orders of magnitude larger 
than all others, and occurs 26 hours before the X flare. The 
next 2 largest spikes occur ~ 38 — 44 hours before the X 
flare. For AR 1261 the 6 largest spikes occur ~ 18 — 39 
hours before the M flare, and all but one are more than 
an order of magnitude larger than the 7th largest spike. 
For AR 1283 the 3 largest spikes occur ~ 22 — 25 hours 
before the X flare, and the 2 largest of these are more than 
an order of magnitude larger than the 4th largest spike. 
For ARs 1429 and 1430, which are magnetically coupled in 
the sense they are merging during the time series, the 9 
largest values of Q occur about one day before the two X 
flares near t=0, but after the X1.1 flare. The meaning of 
the timing of these spikes for ARs 1429/1430, and for AR 
1890 is complicated by the fact that the largest spikes occur 
between X flares. For AR 1890, for values above background 
values, which are © 103 ergs-cm~!-s~!, Q increases from 
the left towards the first X1.1 flare, attains its largest value 
between the two X1.1 flares, and tends to decrease after 
the second X1.1 flare. For AR 2017 the largest spike by an 


order of magnitude occurs ~ 4 hours before the X flare, and 
the next 2 largest spikes occur ~ 90 — 105 hours before the 
X flare. 
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Figure 3. Q for 4 of 7 SF ARs. For AR 1283 the X1.8 and C1.6 
flares overlap. 
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Figure 4. Q for the remaining 3 SF ARs, and N(Q) for all ARs 
(see §V). 


V. THE CDFs OF @ AND FLARES, AND THE 
POSSIBILITY OF SELF ORGANIZED CRITICALITY 
(SOC) IN THE PHOTOSPHERE 
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Figure 5. N(Q) for all C ARs, and 3 SF ARs. 
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Figure 6. N(Q) for remaining 4 SF ARs. 


The CDFs of Q for all ARs, all C ARs, and the 7 individual 
SF ARs are shown in Figs. 4-6. Log-log plots are used 
since if N(Q) = AQ~* where A and s are constant over 
some range of Q, then log N(Q) = —slog@ + logA over 
this range, which is a line with slope —s. CDFs with these 
properties are called scale invariant power law distributions 
since the value of s does not depend on any physical scale 
over this range of Q, so that N(kQ) = constant x N(Q), 
where & is any constant. The lines that are fit to the data 
are generated by the Matlab polyfit function over the range 
of Q indicated in each figure. In most cases the linear fits 
are very good over the indicated ranges. In some cases 
the linear behavior breaks down for N(Q) < 10. Visual 
inspection suggests that in most cases this breakdown is 


due to an insufficient number of data points for the higher 
values of Q, in which case the linear behavior would extend 
to higher values of @ if longer time series were used to 
provide better statistics. However, for SF AR 1166 in Fig. 
5 there is a data point at Q = 1079 ergs-cm~!-s~!, about 
3.5 orders of magnitude to the right of the bottom of the 
linear fit. This suggests that for this AR the scale invariant 


behavior breaks down for Q ~ 1026 ergs-cm7!-s~!. 


The figures, including those of N(Q) for the individual C 
ARs that are not included here, show that the CDF's of the 
14 ARs exhibit scale invariant behavior over ranges of Q 
that extend over ~ 2 — 2.9 orders of magnitude for 3 ARs, 
and over ~ 3.3 — 6.3 orders of magnitude for the remaining 
11 ARs. The CDF's of all ARs, all SF ARs, and all C ARs 
show scale invariant behavior over a range of @ that extends 
over ~ 4—5 orders of magnitude. This behavior is evidence 
that whatever process generates Q above an AR dependent 
threshold value remains the same over the corresponding 
range of Q. This behavior is a necessary but not sufficient 
condition for the system, here an AR NLR, to be in an SOC 
state (Watkins et al. 2016). SOC states arise naturally 
during the evolution of dynamical systems with extended 
spatial degrees of freedom, such as the photosphere and 
corona. They are characterized by long range spatial order, 
are stable in a narrow range of system parameter values, 
and exhibit a high degree of chaos and noise (Bak et al. 
1987; Tang et al. 1987; Kadanoff 1991; Drazin 1992; Bak 
1996; Newman 2005; Watkins et al. 2016), where chaos is 
deterministic but often indistinguishable from noise due to 
insufficient resolution. The transition of a system into an 
SOC state is similar to a phase transition in that the system 
makes a transition from a state without long range order to 
one with long range order. 


For the SF ARs, s =(0.4472, 0.4579, 0.4708, 0.4742, 
0.4858, 0.5101, 0.5143). For the C ARs, s =(0.3351, 0.4080, 
0.4303, 0.4356, 0.4459, 0.5148, 0.5385). The range of s 
is ~ 0.34 — 0.54, with all but one value in the range 
0.4 — 0.54. The mean, median, and standard deviation 
of s for all SF ARs, all C ARs, and all ARs are, respec- 
tively, (0.4800, 0.4742, 0.0252), (0.4440, 0.4356, 0.0675), 
and (0.4620, 0.4643, 0.0524). This shows there is little sta- 
tistical variation in s among the 14 ARs, so s is largely 
independent of individual AR properties. 


The scale invariant behavior of N(Q), and the range, 
mean values, and weak dependence of s on AR are simi- 
lar to the observed properties of the CDF N(£) for coronal 
flares, where F is the total energy released, defined as the 
total amount of magnetic energy converted into particle en- 
ergy. It is observed that above an AR dependent threshold 
of FE, N(E) = kE~° where k is an AR dependent con- 
stant, and ag varies little from one AR to another, being 
largely independent of the differences between ARs such as 
area, total unsigned magnetic flux, and sunspot distribution 
(Datlowe et al. 1974; Wheatland 2000, 2010). Observation 


based estimates of ag have been made for over 40 years 
(Datlow et al. 1974; Hudson 1991; Crosby et al. 1993; 
Shimizu 1995; Aschwanden & Parnell 2002; Aschwanden 
2012, 2013, 2016). For HXR, SXR, and + ray based con- 
structions of N(F), it is found that ag ~ 0.40 — 0.88, with 
almost all values in the range of ~ 0.4 — 0.6. These ob- 
servation based values for ag are consistent with values of 
Qg predicted by first principles, SOC theories, which give 
ap ~ 0.4—0.67 (Aschwanden & Parnell 2002; Aschwanden 
2012, 2013, 2016), and with values of ag predicted by sim- 
ulations of the solutions to lattice based avalanche models, 
which give ag ~ 0.40 — 0.57 (Lu & Hamilton 1991; Lu et 
al. 1993; Charbonneau et al. 2001; McIntosh et al. 2002; 
Aschwanden & Parnell 2002; Aschwanden 2012). 

Therefore the range of s largely overlaps with the range 
of ag, and s and ag both show relatively little variation 
with AR. This suggests a common origin of the photospheric 
N(Q) and the coronal N(£). If they do have a common ori- 
gin, a question is whether SOC theory can help clarify it. 
The first application of SOC theory to propose a physical 
process that gives rise to N(E) is presented in Lu et al. 
(1991, 1993). There it is proposed that the solar corona 
is in an SOC state, and that flares consist of a time series 
of random, spatially distributed avalanche of sub-resolution 
magnetic reconnection events that trigger one another. If 
this theory of flares as a coronal process is correct, a ques- 
tion is whether the photosphere is also in an SOC state, and 
exhibits similar avalanche type flaring events with a simi- 
lar CDF on smaller energy and spatial scales. The form of 
N(Q) found here for 14 ARs suggests the answer to this 
question might be affirmative. The similarity of N(£) and 
N(Q) discussed above does not prove ARs enter an SOC 
state when Q exceeds a threshold value, but it is evidence 
in support of this possibility, and implies it is important to 
further explore this possibility using time series from more 
ARs for better statistics. 


VI. CONCLUSIONS 


The spurious 6, 12, and 24 hour periods in the HMI B 
time series can introduce significant error into B, and quan- 
tities derived from it. 

The CDFs of coronal flares and the photospheric NRL 
resistive heating events predicted by the model presented 
here are very similar in that they are both scale invariant, 
and have power law index ranges that strongly overlap and 
are largely independent of AR. This suggests that the basic 
process that drives coronal flares, and the one that drives 
photospheric heating events are the same process operat- 
ing on two largely different sets of spatial and temporal 
scales, and that this process is part of the evolution of the 
corona and photosphere into SOC states characterized by 
avalanches in the rate of conversion of magnetic energy to 
particle energy. 

The largest photospheric heating events predicted by the 
model, corresponding to the largest spikes in Q for ARs 


with M or X flares are plausibly correlated in time with 
the subsequent occurrence of M or X flares several hours to 
several days later, but the sample size of 14 ARs is too small 
to consider this possibility more than plausible. Analysis of 
time series from more ARs is necessary to determine if such 
a correlation exists, and if it is useful for flare forecasting. 
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